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M  A8SACHUSETTS 


ABSTRACT 


in  this  report  the  detection  and  estimation  of  closely 

.  •.  -  —i  / 

spaced  optical  targets  are  studied  using  simulation.  The 
observed  signal  which  originates  from  two  point  targets  may  ex- 
hibit  only  one  apparent  peak  when  they  are  located  within  one 
detector  width.  The  Akaike  information  criterion  and  a  maximum 
likelihood  estimator  are  used  to  deteat  and  estimate  such  unre¬ 
solved  targets.  For  target  separation  between  3/4  and  1  detector 
width  the  detection  rate  is  high,  the  estimator  is  unbiased  and 
the  estimation  variance  is  close  to  the  Cramer- Rao  bound.  The 
performance  degrades  greatly  when  the  separation  becomes  smaller. 
This  loss  in  performance  is  attributed  to  the  increasing  inter¬ 
ference  between  the  two  targets  and  the  difficulty  in  providing 
a  "good"  initial  guess  for  the  estimator. 
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I. 


INTRODUCTION 


In  the  hierarchy  of  BMD  (Ballistic  Missile  Defense)  systems 
functions,  closely  spaced  object  (CSO)  resolution  oocurs  early  in 
the  sequence  cf  events  and  consequently  influences  the  performance 
of  the  subsequent  functions  to  different  degrees*  It  is  very  im¬ 
portant  to  assess  the  CSO  resolution  capability  of  the  sansor 
systems  employed  in  every  BMD  system  study  before  one  can  determine 
the  overall  BMD  system  performance*  The  CSO  resolution  capability 
is  clearly  dependent  upon  the  sensor  system  and  the  threat  charac¬ 
teristics  considered  in  the  BMD  system  study.  Current  attention  has 
been  focused  on  a  variety  of  passive  optical  sensor  systems  employed 
in  the  Layered  Defense  System  against  threats  at  long  ranges  with 
high  angular  density.  Several  studies  have  been  completed  re¬ 
cently  in  assessing  the  CSO  resolution  performance  for  various 
optloal  systems  [1-4].  The  CSO  parameter  estimation  performance 
can  either  be  predicted  by  theoretical  lower  bounds,  say  the 
Cramer-Rao  lower  bounds  [1-4] ,  or  be  evaluated  by  the  Monte- 
Carlo  simulation  of  specific  algorithms.  In  the  earlier 
studies  [1-4],  the  estimation  accuracy  for  the  intensity  and 
position  of  the  CSO's  were  presented  for  various  lens  apertures 
and  noise  models  under  the  assumption  that  the  exact  number  of 
targets  present  is  known.  The  CSO  detection  performance  was  not 
presented  in  these  studies. 

It  is  the  purpose  of  this  report  to  address  the  following 


1 


issues  through  a  simulation  study  j  is  the  above-mentioned 
theoretical  lower  bound  achievable  in  practioe?  Can  the  number 
of  targets  present  in  the  CSO  cluster  be  determined  with  certainty 
so  that  the  assumption  made  is  true?  In  the  simulation,  a  max¬ 
imum  likelihood  estimator  is  implemented  for  the  CSO  parameter 

I 

estimation  and  the  Akaike  information  criterion  (AIC)  is  em¬ 
ployed  to  determine  the  number  of  targets.  A  specific  sensor  « 

and  noise  model  as  well  as  the  detector  scanning  mode  is  selected 
for  this  study.  For  other  sensors  and  noise  models  and  detector 
patterns,  a  similar  CSO  detection  and  estimation  algorithm  can 
be  implemented  very  easily.  This  work  is  in  its  initial  stages. 

The  findings  reported  here  are  interesting  but  not  necessarily 
complete  and  conclusive.  Further  investigations  are  currently 
in  progress  and  will  be  reported  in  future  reports. 

The  problems  concerned  in  this  study  are  first  stated  in 
section  2.  The  models  of  signal  and  noise  in  a  single  scanning 
detuctor  environment  are  outlined  in  this  section.  The  methods 
for  detection  and  estimation  are  described  in  section  3. 

Section  4  presents  the  Monte-Carlo  simulation  results  of  detection 

and  estimation  performance.  The  estimation  performance  thus  ' 

obtained  is  also  compared  with  the  theoretical  result.  Some 

v 

details  of  the  computational  aspects  of  the  algorithm  and  the 
program  listing  are  attached  in  the  appendices. 
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II.  PROBLEM  STATEMENT 

For  the  purpose  of  this  simulation  study,  the  optical  sensor 
system  can  be  simplified  and  is  represented  in  Fig.  1.  The  opti¬ 
cal  point  targets  which  are  in  the  field-of-view  of  the  sensor 
will  form  an  image  on  the  focal  plane. l, This  image  is  often  re¬ 
ferred  to  as  the  point  spread  function  (PSF) .  An  optioal  detector 
is  usually  employed  to  scan  the  image  in  a  fixed  direction.  The 
spatial  structure  of  the  optioal  image  is  thus  converted  into  a 
temporal  electrical  signal.  This  signal  is  then  subject  to 
amplification,  filtering  and  analog-to-digital  conversion  before 
entering  the  signal  processor.  Noise  sources  which  aan  be  in¬ 
troduced  at  various  points  of  the  system  include  the  background 
radiation  noise,  scanning  noise ,  optical-electrical  conversion 
noise,  amplifier  noise  and  quantization  noise,  For  simplicity, 
it  is  assumed  in  thito  study  that  these  noise  sources  can  be  lumped 
together  and  represented  by  an  additive  white  gausaian  noise  (WON), 
n(t).  The  observations  available  to  the  detection  and  estimation 
processor  can  then  be  written  as 

y(tjj)  “  +  n(tj)  i*l,  2, . . .  ,k  (1) 

where  s^(t)  is  the  desired  signal.  Oiven  these  observations  the 
processor  is  then  required  to  perform  the  following  two  functions i 

1.  Determine  how  many  targets  are  embedded  in  the 
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observed  signal  (detection  problem) 

2.  Obtain  intensities  and  positions  of  the  targets 
(estimation  problem) 

It  is  the  purpose  of  this  report  to  describe  an  algorithm  for 
such  a  processor  and  to  evaluate  its  performance. 

Suppose  that  there  are  n  point  targets  which  lie  along  the 
scanning  direction  of  the  detector. *  The  desired  signal  sd(t) 
is  simply  given  by 

n, 

ad(t)  •  ^  (2) 

where  a^  and  are  the  intensity  and  position  of  the  ith  target; 
aQ(t)  is  the  basic  target  response  generated  by  a  target  of  unit 
intensity  lying  on  the  optical  axis.  The  shape  of  s0(t)  depends 
upon  the  FSF  of  the  particular  aperture  and  the  scanning  response 
function  of  the  detector.  Suppose  the  aperture  is  annular  with 
50t  obscuration;  the  detector  response  function  is  uniform  and 
equal  to  unity  within  a  rectangular  gate  and  equal  to  zero 
elsewhere.  Then  sQ(t)  is  given  by  [6). 

>.fly/2  -at+9x/2 

s0(t)  “  /  /  ef(x,y)dxdy  (fla) 

•/-0y/2*/  at-flx/2 

"This  assumption  lias  to  be  mads  beoause  the  cross-scan  position 
of  a  target  can  not  be  resolved  by  a  single  detector. 
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where 


2 


8  J  [Tr(x2+y2))l]  4  J1[5-(x2+y2)lS]i 

af(x,y)  -  -  — — H - LA.  ......  ......  ,  (3b) 

3  ii  (x^+y2 ) ~  3  1r  (y2+y2)  »  j 


Hare,  is  the  Bessel  function  of  first  kind;  a  is  the  an- 

guiar  scanning  rate  normalized  by  the  optical  diffraction  limit 
\/&t  8X  and  are  the  in-scan  and  croes-soan  angular  dimensions 
of  the  detector  normalized  by  1/d.  Without  loss  of  generality, 
the  scanning  rate  oan  be  assumed  equal  to  1  and  thus  the  time 
variable  t  is  equivalent  to  the  angular  variable  9.  These  two 
variables  will  be  used  interchangeably  throughout  this  report. 

The  basic  optical  pulse,  s0(t),  is  depicted  in  Fig.  2  for  a 
detector  with  8X"2  and  8y«6.  Note  that  there  is  a  slight  overshoot 
near  the  center  of  the  pulse.  Fig.  3  illustrates  several  sample 
waveforms,  y(t),  from  detectors  of  identical  size  for  the  aase  of 
2  CSO'a  with  various  intensities*  and  positions .  Tho  variance  of 
the  WGN  is  equal  to  1,  it  can  be  seen  that  the  2  targets  which  are 
separated  by  2.5  X/d  (Fig,  3a)  correspond  to  two-peaks  of  y(t). 

In  this  case  they  can  be  detected  and  estimated  with  a  matched 
filter  followed  by  a  peak  detector.  However,  in  Figs.  3{b)  -  3(d), 
the  two  targets  which  are  located  within  one  detector  width (2X/d) 
interfere  with  each  other  to  such  a  degree  that  only  on^  apparent 


as  the  peak  at  the  center 


to  10  for  a  *  9, OB. 
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Fig.  3(a-d).  Some  examples  of  the  observed  signal,  y(t),  in  the 
case  of  two  CSO's>  the  variance  of  the  WON  ia  equal  to  1. 
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peak  is  observed  in  the  resulting  noisy  waveform.  A  simple  peak 
detector  may  not  be  capable  of  resolving  satisfactorily  the  two 
targets  in  these  examples.  Other  more  sophisticated  detection/ 
estimation  schemes  are  required  for  this  purpose. 
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Ill.  METHODS 


3.1  Akaike  Information  Criterion 

One  approach  for  determining  the  number  of  target*  imbedded 
in  the  observed  data  is  to  apply  the  generalised  likelihood  ratio 
(GLR)  test.  There  are  some  difficulties  in  applying  this  method. 
First,  the  distribution  of  the  GLR  is  hard  to  find  so  that  the 
behavior  of  the  test  may  not  be  known  exactly.  Second,  since 
the  test  can  only  be  applied  to  two  desses  at  one  time,  multiple 
application  of  the  test  is  required  for  the  present  problem. 

Third,  the  choice  of  the  test  threshold  is  usually  very  subjective. 

Akaike  has  advooated  a  new  approach,  termed  the  Akaike  in¬ 
formation  criterion  (AIC) ,  for  determining  the  order  of  the  correct 
model  (for  the  problem  here,  the  order  is  the  number  of  targets 
present)  [5].  This  information  criterion  is  based  on  an  extension 
of  the  maximum  likelihood  principle  starting  from  the  fundamental 
notation  of  entropy  in  statistical  mechanics  and  the  Kullbaok- 
Leibler  information  quantity.  The  final  statistic  used  to  optimally 
choose  the  order  is  defined  by 


Aic  «  (-2)  log  (maximum  likelihood)  +  2 (number  of  free  parameters) 
0 

(4) 

The  correct  model  is  that  which  minimises  this  criterion. 
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Note  that  the  first  term  in  the  definition  of  the  AX C  represents 
a  penalty  of  "poor  fit"  and  the  second  term  characterizes  in¬ 
creased  unreliability.  This  second  term  is  essential  because  the 
maximum  value  of  the  likelihood  function  with  the  higher  order 
model  (model  with  more  parameters)  is  usually  greater  than  that  v 
of  the  smaller  model  and  therefore  without  this  term/  the  model 
with  higher  order  would  be  favored.  Qualitatively  speaking,  the 
AIC  provides  a  mathematical  formulation  of  the  principle  of 
parsimony  in  model  building. 

The  AIC  has  found  many  applications  in  various  fields,  par¬ 
ticularly  in  the  autoregressive  model  fitting  of  time-series 
analysis.  Sinoe  the  theory  is  general,  it  is  easily  applied  to 
the  detection  problem  stated  in  the  last  section,  Two  obvious 
advantages  for  using  the  AXC  in  this  problem  can  be  seen  from 
Eq.  (4) t 

1,  The  AIC  is  easy  to  apply;  Eq,  (4)  is  very  simple. 

2.  The  AXC  combines  the  detection  problem  with  the  estima¬ 
tion  problem.  For  detection  (determining  the  order  of 
the  model)  the  information  from  estimation  (finding  the 
maximum  likelihood  function)  is  needed.  Once  the  de- 
tection  is  done,  the  maximum  likelihood  estimates  of 
parameters  in  the  correct  model  are  also  available  with1 
out  additional  effort. 

Although  the  AXC  is  derived  from  the  ideas  of  information 


theory,  there  seems  to  be  ho  particular  basis  for  the  penalty 
factor  2.  Some  Investigators,  applying  the  AIC  in  their  Par¬ 
ticular  problems,  reported  that  this  factor  should  be  between 
3.5  and  4  [7],  it  can  be  shown  that  the  AXC  decision  rule  is 
equivalent  to  the  hypothesis  testing  procedure  at  an  appropriate 
significance  level  In  the  bails'  of  ''two  classes. :  ,tising  different 
values  for  the  penalty  factor  is  analogous  to  adjusting  the 
significance  level.  The  effect  of  this  factor  will  be  reported 
in  this  study.  For  this  purpose,  the  AXC  oan  be  rewritten  as 

AIC(i)  «  (-2)  lc<7e (maximum  likelilioud)  -i  'i(2i)  (5) 

where  n  is  the  penalty  factor  and  i  is  the  number  of  targets  pre¬ 


sent. 


Xt  is  reasonable  to  assume  that  AXC(i)  is  a  discrete  uni- 


modal  funotion  of  i  for  a  fixed  value  of  n.  Therefore  the  defec¬ 
tion  procedure  is  simply  as  follows t 

Step  li  Start  with  i«0>  compute  AXC  ■  AXC(O). 

Step  2t  Increment  i  by  1. 

Step  3t  Compute  AXC(i). 

Step  4 1  Xf  AXC (i)  is  greater  than  AXC,  go  to  Step  5.  Other* 
wise,  set  AXC  «  AXC (i)  and  go  back  to  Step  2. 


target. 


position,  are  associated  with  each 
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Step  5 i  Stop;  the  number  of  targets  present  is  equal  to 
(i-1). 

.  ..  .  -  {>  •  .  ’  '  '  ’ 

3.2  Maximum  Likelihood  Estimator 

.The  maximum  likelihood  estimator  (MtE)  la  implied  in  the  AIC 
procedure.  This  estimator  posesses  several  nice  properties.  It 
can  be  shown  that  the  MLS,  under  rather  general  conditions,  is 
asymptotically  unbiased  And  efficient;  it  yields  the  same  results 
as  the  least-square  method  for  the  case  of  additive  white 
gaussian  holes. 

For  the  signal  modal  described  in  section  2,  the  likelihood 
function  oan  be  written  as 


where 


&  (x ,  c) 


irrery  •**(- <ro„a>  Vv.>) 


(2w>-  .  c 


(6) 


°n  "I 


/vW*  •  •  •  •  *Bo(tr'rn) 


(7) 


’WV- 


■o 


w  (V'V  Tl,,,'l'n)  xn+l"*,:<2n) 


x 


J(x,o)  ■  log#4(x,a)  « 

-  -  y  loge(27r)-klogao  -  j^r<£*Qna)T(£-Qna)  <*) 

The  maximum  value  of  this  function  i*  the  first  term  needed  to 
evaluate  the  AIC. 

It  ia  uaually  not  naaaaaary  to  estimate  tha  unknown  a  ex¬ 
plicitly  and  it  can  be  dropped  from  tha  expreaeion  for  J.  By 
taking  tha  derivative  of  J  with  reapeot  to  a  and  aetting  it  to 
aero,  9  ia  obtained  aa 


By  replacing  o  with  a,  Bq.  (8)  becomes 

J(x)  -  {logoff)  +  1  +  l°Se(jjr(£-an£)T(£-Qn*)  (10) 

Using  this  logarithmic  likelihood  function  with  tha  firat  two 
nuisanoa  terms  discarded,  the  AIC  given  by  Eqt  (5)  can  ba 
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written  as 


AIC(i)  *  klog.c2*  2ni.  (11) 

© 

It  is  clear  that  as  far  as  the  maximisation  is  oonoerned 
J(x)  in  Eq.  (10)  can  be  replaced  with 

J(x)  -  “ (£-Qn£) T (X-Qn*) •  (12) 

By  setting  the  gradient  of  J(&)  with  respect  to  ft  equal  to  zero, 
the  intenatiy  estimate,  a,  can  be  obtained  as 

a  -  (on\)  an\  •  (i3) 

substituting  ft  in  Eq.  (12),  J  oan  be  rewritten  as 

J<I>  -  lCT[Qn(0nT0n)  lOnT-l]i:  .  <1<! 

The  maximization  of  the  likelihood  function  oan  be  done  with  res¬ 
pect  to  expressions  given  by  either  Eq,  (12)  or  Eq,  (14),  The 
former  involves  2n  parameters  and  the  latter  involves  only  n 
parameters  but  requries  ratrix  compuatione  which  usually  need 
more  computer  time.  Experience  seems  to  indicate  that  it  is  easier 
to  use  Eq,  (12) , 
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It  should  be  noted  that  the  J  funotion  Is  nonlinear  on  the  un¬ 
known  parameters.  Its  maximization  is  implemented  here  be  using 
the  Quasi-Newton  method  [8]  which  is  an  iteration  procedure  start¬ 
ing  with  an  initial  guess.  This  method  is  appealing  because  it 
possesses  the  quadratic  convergence  property  near  the  maximum  of 
the  criterion  fuhotion  (like  the  Newton  method)  and  avoids  the 


difficulty  involved  in  computing  the  inverse  of  the  Hessian 

matrix  at  each  iteration  (unlike  the  Newton  method) .  At  eaah 

3J(*) 

iteration,  the  gradient,  VJ(x)  ■  {  ■  .  .  i-l,.,.,2n),  is  required, 

3xi 

which,  from  Eq.  (12),  is  given  by 


Here,  from  Bq.  (2) , 


t 


* 


3sj  (t)  80(t-*fj)  X,»Sj 

— 2 -  “  ,  (16) 

3xi  -*i*o^-Ti>  Xi*’ri 

and,  from  Eq,  (3) , 

&Q(t)  “  a  f ^  [sf  (c*t+3x/2,  y)  -  sfi(ctt-$x/2,  y)jdy.  (17) 

The  unknown  parameters  are  not  entirely  free  but  subject  to 
two  types  of  constraints.  The  intensity  of  any  target,  a^,  is 

physically  constrained  to  be  non-negative.  This  constraints  can 
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2 

be  released  by  substituting  b?  for  a^.  By  doing  so,  nothing  is 
changed  except  that  the  parameter  sat  becomes  x  -  (b^,.. ,,bn> 
T1f.*.#Tn}T  and  Eq.  (16)  is  replaced  by 


The  final  estimate  a^  can  be  obtained  by  setting  a^  ■  £;2. 

Suppose  now  the  sequenoe  of  observations,  £,  is  obtained  in 
the  range  from  90-y/2  to  9q+y/2  where  y  is  the  angular  distance 
covered  by  the  observations  and  8Q  is  the  central  point  of  the 
range.  In  preotiae,  it  can  be  assumed  that  the  admissible  range 
for  t^,  center  of  the  target  response,  is  also  within  this 
range,  in  other  words,  only  targets  with  center  falling  in1 this,  . 
range  will  be  identified.  This  additional  range  constraint  will 
be  incorporated  into  the  Quasi-Newton  method. 

The  choice  of  initial  guess  is  a  crucial  step  for  the  Quaai- 
Newton  method,  A  good  initial  guess  can  lead  the  iteration  process 
to  converging  to  the  correct  maximum  in  a  relatively  small  number  of 
iterations.  On  the  other  hend,  for  a  bad  initial  guess,  the  itera¬ 
tion  process  might  converge  to  a.  local  maximum,  oscillate  around 
the  maximum  or  not  converge  at  all. 

One  method  for  generating  the  initial  guess  is  the  pure  ran- 
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* 

dom  search.  This  method  consists  of  computing  J(,£)  at  N  random 
points  drawn  from  a  probability  distribution  uniform  over  the 
entire  parameter  space  and  selecting  the  point  with,  the  greatest 
value  of  J(x)  as  the  initial  guess,  if  it  is  assumed  that  eaah 
parameter  oan  vary  between  0  and  1001  and. that  the  optimal  para¬ 
meter  set  which  correspond  to  the  global  maximum  of  J(V)  is  to  be 
located  within  101  of  each  parameter  than  the  probability  of 
looating  the  optimum  in  N  trials  is  [9] 

N 

p  -  l-(l-10“n)  .  (19) 

Conversely  the  number  of  trials  required  to  have  a  90*  probability 
of  locating  the  optimum  is 

N  -  1/log (—i—)  *  2.3xlOn.  (20) 

■  iriO “n  _■  :  •; ..  7  : 

Obviously,  the  required  number  of  trials  increases  rapidly  with 

the  number  of  unknown  parameters,  n.  if  there  exists  a  single 

target  (n»l)  the  pure  random  search1  Seems  able  to  yield  a  "good" 

initial  guess  within  a  reasonable, computation  time.;  However 

for  more  than  one  target  this  method  becomes,  impractical. 

A  more  practical  approach  is  to  use  the  pure  random  search 

in  conjunction  with  a  priori  knowledge*  As  pointed  put  earlier, 

*!For~tHe~  purpose  of  select ing  an  initial  guess,  Bq.  (14)  instasd 
of  Bq.  (12)  is  used. 


IB 


in  applying  the  AIC-  procedure  to  determine  the  number  of  targets , 
the  statistics  Aic(i),  i-'l,  2,...,are  computed  in  ascending 
sequence.  At  each  step  the. likelihood  function  is  maximized  iter¬ 
atively  starting  with  an  initial  guess.  It  seems  feasible  to 
select  the.  initial.  g,uesd\o£  'the;.  (i+l)th  step  -in  suoh  a  way  that 
the  initial  values  of  the  first  i  parameters  are  equal  to  the  i 
estimates  from  the  previous  step  and  the  initial  value  of  the 
remaining  parameter  is  ohosen  from  pure  random  searoh.  Thus,  the 
initial  guess  for  any  model  (number  of  target)  is  obtainod  by  per¬ 
forming  a  one-parameter  searoh  whioh  is  easier  numerically. 

This  approach  guarantees  that  the  maximum  likelihood  function 
of  the  (i+l)th  step  is  no  less  than  that  at  the  ith  step. 

However,  if,  for  the  model  of  i+1  target*,  the  first  i  components 
of  the  true  optimum  in  the  i+1  parameter  spaas  ere  far  from  the 
optimal  point  ..determined  in  step  i,  the  initial  guess  obtained 

W 

in  this  way  is  usually  not  "good"  for  the  (i+1)  Step,  The 
actual  implementation  of  this  approach  is  described  in  the 
appendices. 
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IV,  SIMULATION  -  RESUI/CB.  AND  DISCUSSIONS 

The  performance  of  detecting  the  number  of  targets  present  and 

estimating  the  parameters  of  these  targets  is  evaluated  by  Monte- 

■  l,« 

Carlo  simulation.,  The  simulation  algorithm  has  been  described  in 
the  previous  section'  and  :more  details  are  given  in  the  appendices'. 

In  presenting  the  results,  some  system  parameters  must  be  specified. 
It  is  assumed  that  the  signal  available  to  the  processor  is  observed 
in  an  angular  range  of  +  6,3  A/d  from  the  center  of  the  focal  plane. 
The  signal  is  sampled  uniformly  in  this  range  at  an  interval  of 
.2  A/d.  .  The  total  number  of  observations  is  64.  The  signal  is 
the  product  of  an  annular  aperture  with  50%  obsauration  and  a 
scanning  photo-detector  with  in-scan  and  oross-scan  dimensions 
equal  to  2  A/d*  and  6  A/d  respectively.  The  number  of  Monte-Carlo 
runs  is  fixed  at  100  for  every  cane  discussed  in  the  following. 

4,1  Detection  Performance 

First  consider  the  case  where  at  most  one  target  may  axiet. 

One  must  decide  between  two  hypotheses!  Hqi  no-target  and  H^i 
one-target  present.  For  this  case  it  is  only  necessary  to  compute 
AXC(O)  and  AXC (1) ,  and  the  decision  rule  accepts  H0(H^)  if  AXC(O) 
(aic (1) )  is  the  smaller  of  the  two.  The  detection  performance 

•tfhe  in-scan  dimension ’"of" the  detector  is  approximately  equal  to 
the  diameter  of  the  blur  size  (diameter  of  the  first  dark  ring)  of 
the  PSF, 
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curve  (false  alarm  rate  vs.  leakage  rate)  ie  shown  In  Fig,  4. 

The  false  alarm  rate  (P£)  is  the  probability  of  aocepting  given 
that  Hq  ie  true  and  the  leakage  rate  (P^)  is  the  probability  of 
accepting  H0  given  that  is  true.  The  parameter  d  in  the  figure 
is  signal-to-noise  ratio  related  and  defined  as 


Here  a  is  set  equal  to  1  and  a  varies  among  .908,  .454.  .227  and 
.114  for  the  four  curves  shown.  Bach  point  on  the  performance 
curve  corresponds  to  one  value  of  the  penalty  faotor  n  in  the 
A1C  defined  by  Bq.  (5).  The  value  of  n  is  established  by  the 
aatual  operating  point  which  is  determined  by  the  requirements 
for  and  P£, 

Now  consider  the  case  where  there  are  two  CSO's  (n»2)  sepa¬ 
rated  by  A9  which  is  less  than  or  equal  to  2X/d  (one  detector 
width).  It  is  always  assumed  that  the  two  targets  are  located  at 
TjU-Ae/2  and  t2  “A0/2 ,  and  o  is  equal  to  1.  Fig.  5  shows  the  pro¬ 
babilities  of  identifying  correctly,  from  a  given  observation,  2 
targets,  P(2),  less  than  2  targets,  P(<2),  and  more  than  2  targets, 
P(>2).  The  observation  comes  from  two  equally  strong  targets  with 
signal-to-noise  ratio  (SNR  )  of  10.  The  penalty  factor,  n  is 
chosen  as  3  hare.  The  curves  shown  arc  drawn  by  connecting  the 


"The  “SUIT  of  a  target  is  defined  as  the  intensity  of  the  target  at 
the  center  of  its  pulse  shape  divided  by  the  RMS  of  the  noise. 

The  interference  noise  is  not  included. 
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T»rg«t  g»p*rition 
In  Unit*  of  Vd 

Fig.  S.  Probabilities  of  identifying  2  targets  P(2),  less 
than  2  targets  P(<2),  and  more  than  2  targets  P(>2),  while 
the  separation  between  the  actual  targets  is  varied. 


finite  simulation  points.  The  probability  of  oorreotly  identifying 
2  targets  is  as  high  as  90%  for  separation,  Ad,  of  2X/d  (one 
detector  width).  However  it  is  less  than  25%  for  separations 
smaller  than  ,5X/d  (quarter  of  the  da  tea  tor  width).  The  error 
is  mainly  because  the  2  CSO's  are  identified  as  a  single  target. 

For  separation,  A6,  between  *5X/d  and  X/d,  the  correct  detection 
rate  is  between  25%  and  65%,  and  mis-datection  is  attributed  to 
both  under  and  over  identifying  the  number  of  targets.  The  data 
at  no  reparation  (A 0-0)  can  be  viewed  from  another  angle.  At 
A 9-0,  the  2  CSO's  are  coincident  and  indistinguishable  from  a 
single  target.  Therefore,  the  data  indicates  that  a  single  target 
with  SNR-20  has  95%  probability  of  being  identified  as  a  single 
target  and  5%  probability  as  2  CSO's. 

Fig,  6  illustrates  the  effect  of  choosing  different  values 
for  the  penalty  factor  n.  It  can  be  seen  that  using  larger  values 
of  n  can  inarease  the  detection  rate  for  separation  A 6> X/d,  but, 
at  the  expense  of  poorer  performance  for  A9<X/d,  There  appears  to 
be  no  obvious  way  to  select  n  optimally.  This  does  not  oontradict 
Akaike'a  theory  because  no  proof  for  the  optimality  of  using  2 
for  n,  as  proposed  by  Akaike,  has  been  given.  Although  not  shown 
in  the  figure,  it  is  worthwile  to  point  out  that  for  n"0,  P(2)  is 
equal  to  0%  for  any  separation,  A6.  In  other  words,  the  number 
of  targets  present  is  never  correctly  identified  as  2  when  the 
penalty  term  of  the  A1C  statistics  vanishes. 
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T«*g*t  Itpuition 
in  Unit*  of  k/d  otr  dttwtor  width 


f'iflf.  6 .  Effect  of  ehooaing 
penalty  factor  in  Bq,  (11), 


diffarant  valuaa  for  the 
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Fig.  7  shows  ?(2)  as  a  function  of  A0  for  several  different 
sets  of  target  intensity.  The  variance  of  noise  remains  equal  to 
1.  it  is  observed  that  the  oorreqt .  detection  rate  is  gainer  ally 
higher  for  greater  intensity  except  for  the  sudden  drop  of  the  top 
curve  at  ,&e*X/d..,  This:4r9p.'in  detection  ret* -i*- due  to  the  feat  . 
that' in  many  of  the  100  Monte  Carlo  rune  for  this  case,  the  ini¬ 
tial  guess  provided  by  the  simulation  algorithms  causes  the 
Quasi-Newton  method  to  oscillate  or  to  converge  to  a  local  maximum 
for  the  2-target  model  but  leads  to  a  final  estimation  for  the  3- 
targat-model,  which  includes  2  targets  close  to  the  true  ones  and 
a  third  target  of  insignifioant  intensity.  Note  that  an  axperimen-  . 
tor-supplied  "good"  initial  guess  for  the  2-target  model  can  bring 
up  the  detection  rate  at  this  point.  As  to  why  the  algorithm 
fails  at  this  particular  point,  no  ,f.atia  factory  explanation  has 
been  found,.;  •  - 

In  Fig.  5  the  true  model  has  2-target  with  a.]«a2-».08  and 
Ti«,-t2"-A6/2.  Suppose  the  initial  guess  of  the  Quasi-Newton  pro¬ 
cedure  is  supplied  by  the  experimentor  instead  of  the  algorithm 
itself  and  the  experimentor  intelligently  selects  s^i*18.16,t^*0 
fcr  th.  initial  «u...  ol  th.  l-t.rj.t  modal, 

T1°»-T2°»-h9/2  for  the  2-target  modal  and.  a^^aj^S.OS,  a3  -.1, 
t1°—T2°«-ao/2,  t^-0  for  the  3-target  model.  The  resulting,  de¬ 
tection  performance  ie  shown  in  Fig.  a,  A  comparison  of  Figs. 

5  and  e  shows  that  tha  "good"  initial  guess  increases  the  correct 
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good"  initial  guaaa 


0  .33  .3  ,78  1.  aK 

Target  Separation 
in  Unit#  of  \/A  or  dateotor  width 

Fig.  8.  Comparison  of  detaction  performance  between 
using  a  "good"  initial  guess  and  using  the  initial 
guess  provided  by  the  simulation  algorithm. 


detection  rate  significantly.  For  separations  as  small  as  .8A/d 
(two  fifths  of  the  detector  width)  the  detection  performance  is 
almost  perfeot.  The  incorrect  detection  for  even  smaller  separa¬ 
tions  is  because  only  one  target  is  reaognized.  The  chance  of 
identifying  more  than  2  targets  is  almost  nil  in  the  entire  range 
of  separations  considered. 

4.2  Estimation  Performance 


The  performance  of  an  estimator  is  usually  evaluated  in  terms 
of  estimation  bias  and  variance.  From  the  N  simulation  runs,  the 
variance  of  an  estimate  is  computed  as 


•  g  *ct) 

.<i  n  jml  i  i 


where  x^  is  the  estimate  of  the  parameter  x^  in  the  run  and 


(23) 


is  the  mean  of  the  estimate.  The  bias  of  the  estimate  is  then 
given  by 


bxt  ■  Xjl  -  xi  .  (24) 

This  sample  variance  can  be  compared  with  the  Cramer- Rao 
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bound  (CRB)  which  Is  a  lower  bound  on  the  varianoe  (or  oovariance 
matrix)  of  an  unbiased  estimator.  In  the  oaae  of  a  tingle  target 
in  additive  white  gauctian  nolle,  the  CRB  for  the  ettimation  of 
parameter  vector  x  can  be  obtained  by 


CfXj.)  -  (P)jJ 


(25) 


aad(V  ^d1^1 
3x^ 


(  26) 


2 

where  F  ii  the  Fiaher  information  matrix,  cr  the  varianoe  of  the 
noita  and  sd(t)  the  deaired  signal.  A  aat  of  diaorete  obaervationa 
is  aaaumed.  Thia  bound  ia  uaually  tight  whan  the  signal-to-noiae 
ratio  la  high. 

Figs.  9  (a)  and  (b)  compare  the  square  roots  of  the  o(x^) 

and  a?  for  the  intenaity  and  portion  respectively,  in  the  oaae 
*L 

of  two  targets  with  SNR' a  both  equal  to  10.  The  angular  separa¬ 
tion  between  the  two  targets  ia  the  control  variable.  The 
simulation  ia  only  run  at  certain  values  of  target  separation. 

For  the  data  shown  in  these  figures  the  deteotion  procedure  ia 
omitted  by  assuming  that  the  number  of  targets  present  is  known 
exactly  in  advance.  That  ia  to  say  a  2-target  model  ia  always 
assumed.  The  initial  guess  for  this  model  la  provided  in  two 
different  ways  for  comparison.  The  first  uses  the  procedure 
described  in  section  3.2  and  the  second  relies  on  an  "intelligent" 
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•*  »•  i-»  a.  t,t  x/t, 

Mptrttion 

in  Unit*  of  dttootor  width  or  \/i 

Fig.  9(a-b),  Estimation  performance  and  comparison* 
(a)  intensity  estimation,  (b)  position  estimation* 
one  example. 
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experimentor  who  select*  the  initial  value  of  the  parameter  in  the 
neighborhood  of  the  true  one.  Mots  that  since  a^aj  and 
|  "t  1 1  ™  |  t  2  i  r  the  concerned  statistics  for  the.  first  target  should 
be  very  close  to  those  for  the  second  target  and  so  no  distinction 
between  them  is  made  in  this  figure..',  .  '  ' 

From  Fig.  9,  it  can  be  seen  that  when  the  initial  guess  is 
provided  by  the  algorithm  automatically,  the  estimate  variance 
agrees  with  the  CRB  very  well  at  target  separations  of  l.SX/d, 

2X/J  and  2.5X/d.  Recall  that  the  observed  signal  exhibits  two 
peaks  for  separation  of  2.5X/d  as  shown  in  Fig.  3(a).  The 
sample  variance  is  significantly  larger  than  the  CRB  for  separa¬ 
tions  less  than  1.2X/d«  When  the  "good"  initial  gueBses  are 
employed,  the  estimation  variance  becomes  quite  close  to  the  CRB 
except  at  the  point  where  separation,  &9«.5X/d,  This  once  again 
demonstrates  the  importance  of  the  initial  guees  in  an  iterative 
optimization  algorithm. 

Figs,  10(a)  and  (b)  show  the  same  types  of  comparison  as 
that  of  Figs.  9(a)  and  (b)  for  two  targets  with  SNR's  both 
equal  to  13,  Similar  behaviors  ars  observed.  It  should  be  pointed 
out  the  yc"(a)/a  and  y  C  ( t )  of  this  example  are  smaller  than  those 
of  the  former  example  by  a  factor  of  1,3,  which  is  exactly  equal 
to  the  ratio  of  the  two  corresponding  SNR's. 

¥In  ^he  'first “case#*  We'1' estimation  procedure  starts  with  the  one- 
targst  model  and  ends  with  the  two-target  model.  In  the  second 
case,  it  is  only  applied  to  the  two-target  model. 
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Based  on  the  same  simulation  as  used  in  Fig.  9,  the  biases 
and  the  RMS  errors*  of  the  intensity  and  position  estimations  are 
shown  in  Fig.  11(a)  and  (b) ,  Here  the  initial  guess  is  generated 
by  the  estimator  itself,  it  can  be  seen  that  the  estimator  is 
actually  biased  in,  the  region  where  the  sample  variance  deviates 

*  . .  .."  i  •  M.  ,•  i  ...  •  i  .  r  ' .  .  -..i 

significantly  from  the  Cramer-Rao  bound.  Meanwhile,  the  intensity 
bias,,  b  ,  is  relatively  small  compared  to  the  corresponding  sample 
standard  deviation  (roughly  by  ah  order)  over  the  entire  CSO 
region.  However,  at  vary  small  separation  (<lA/d) ,  the  position 
bias  is  significant  and  contributive  to  total  RMS  error.  If 
only  those  runs  in  which  the  number  of  targets  is  identified  as 
2  are  used  in  computing  the  statistics,  tha  bias,  sample  varianoa 
and  RMS  error  aan  be  recued  slightly. 


-  * . . —  ■  3 — 

The  mean-square  arror  e  is  tha  sum  of  the  varianoe  and  the 
2  2  2 

squared  bias  e  «o*+b  , 
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4.3  Diacusalon 


In  this  simulation  study/  the  pulse  shape  of  the  response  of 
the  photo-detector  to  a  unit  ■’intensity  point  target  is  assumed 
known  exactly  in  advance.  What  the  detection/estimation  processor 

a ..  '  v, 

does  is  to  adjust  the  intensities  and  positions. of  a  certain 
number  of  ideal  pulse  shapes  to  best  match  this  "synthesized" 
signal  with  the  given  observation.  If  the  pulse  shape  is  dif¬ 
ferent  from  the  one  assumed  in  this  report  due  to  differences  in 
aperture  shapes/  detector  responses/  deteotor  configurations 
eta,,  the  simulation  program  is  still  applicable  as  long  as  the 
pulse  shape  and  its  derivative  can  be  made  available.  However,  the 
performance  of  the  detection  and  estimation  process  may  vary 
significantly  with  different  pulse  shapes, 

The  algorithm  described  in  this  report  is  designed  to  process 
the  CSO  segments  rather  than  the  entire  observation.  It  seems 
feasible  to  use  this  algorithm  as  the  second  stage  of  a  two-stage 
signal  processor  which  first  process  the  original  signal  to  identi¬ 
fy  the  isolated  "resolved"  targets  and  to  separate  them  from  the 
potential  CSO's.  This  approach  may  be  quite  efficient  computa¬ 
tionally  if  the  probability  of  occurrence  for  CSO's  is  much  lower 
than  that  of  isolated  targets. 

One  version  of  the  Quasi-Newton  method  is  implemented  here  for 
the  optimization  of  the  likelihood  function.  There  exist  some 
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other  methods  for  this  purpose ,  which  can  be  divided  into  two 
general  categories;  the  random  method  and  the  gradient  method. 
Usually ,  no  single  method  is  best  for  all  types  of  nonlinear 
optimization  problems.  No  evaluation  of  different  methods  on  the 
present  CSO  problem  is  undertaken  in  this  report.  However  it  is 
the  author's  feeling  that  selecting  a  particular  optimization 
method  may  not  be  as  important  aa  devising  an  efficient  and  premia 
ing  way  to  provide  the  initial  guess. 

A  single  scanning  detector  is  assumed  in  this  report.  Since 
this  detector  can  not  resolve  the  croee-scen  component  of  a  tar¬ 
get's  position,  the  targets  are  assumed  to  lie  along  the  in-scan 
direction.  To  determine  the  in-eoan  and  oross-scan  positions 
of  a  target,  other  detector  configurations  such  as  the  chevron 
(a  pair  of  detectors  oriented  in  different  directions)  should  be 
used. 

In  this  simulation  study,  an  additive  white  guassian  noise  is 
assumed.  This  assumption  is  valid  when  the  thermal  noise  is  the 
dominant  noise  source  in  the  optical  sensor  system.  However  this 
noise  model  becomes  inaccurate  in  the  so-called  shot-noise  limited 
case  where  the  noise  level  is  dependent  on  the  signal  [3],  Sven 
for  this  case,  the  detection/estimation  scheme  presented  in  this 
report  can  apply  except  that  the  likelihood  function  and  its 
gradient  should  be  reformulated.  The  mathematics  involved  is, 
of  course,  more  complicated  but  still  tractable. 


37 


V. 


CONCLUSION 


In  this  report  the  detection  and  estimation  problems  of 
closely  spaaed  optical  point  targets  are  studied  ucing  simulation. 
An  annular  aperture  of  50%  obscuration  and  a  scanning  photo^detec- 
tor  with  in* sc an  and  crotti-scan  dimensions  equal  to  '2X/d  and  6X/d, 
respectively,  are  employed.  The  observed  signal  which  originates 
from  two  point  targets  exhibits  a  single  apparent  peak  when  they 
are  separated  by  less  than  one  deteotor  width.  To  deteot  and 
estimate  suoh  "unresolved"  targets,  the  Akaike  information 
criterion  and  a  particular  maximum  likelihood  estimator  are 
utilized.  The  estimator  is  in  faot  a  nonlinaar  estimation 
algorithm. 

Several  examples  have  been  examined.  It  is  found  that  for 
target  separation  between  3/4  and  1  detector  width,  the  correct 
detection  rata  is  fairly  high,  the  estimator  unbiased  and  the 
sampled  variance  close  to  the  cramer-JRao  bound.  However,  the 
detection  and  estimation  performances  degrade  significantly  for 
smaller  separations.  This  is  unavoidable  because  when  the  two 
targets  get  closer  they  become  more  indistinguishable  from  a  single 
target,  particularly  in  the  presence  of  noise.  More  importantly, 
the  difficulty  in  providing  algorithmically  a  "good"  initial  guess 
for  the  estimator  contributes  to  the  poor  performance.  Un¬ 
doubtedly,  the  performance  of  the  algorithm  can  be  improved  if  a 
more  intelligent  way  of  choosing  the  initial  guess  can  bo  devised. 
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APPENDIX  At  COMPUTATIONAL  ASPECTS 


Programs  have  been  written  in  accordance  with  the  principles 
described  in  section  3  in  .order  to  simulate  the  detection/estimation 
processor.  The  performance  of  the  processor  is  oloseiy  related  to 
the  particular  implementations  adopted,  Therefore,  before  present¬ 
ing  the  simulation  results,  the  program  actually  used  has  to  be 
specified  more  carefully.  The  top-level  flow  charts  and  details 
of  some  subroutines  are  given  in  this  the  following  appendices. 


For  practical  reasons,  the  simulation  is  done  in  two  steps 
by  two  different  main  programs.  The  first  one  (TABLE)  is  used 
to  create  the  standard  pulse  shapes  of  sQ(t)  and  sQ(t)  in  advance 
for  use  in  the  second  one  (CSOMCS)  which  is  the  program  performing 
detection  and  estimation. 

The  top-level  flow  aharts  of  these  two.  programs  are  depicted 

i  ■ 

in  Figs.  Ai(a)  and  (b)  respectively  while  their  listings  are  giver 
in  the  appendix.  In  these  eharts  the  functions  ot  some" blocks 
are  implemented  by  subroutines.  The  names  of  these  subroutines  ' 
are  put  down  beside  the  associated  blocks.  Some  of  them  will  be 
further  explained  in  the  following  subsections.  Note  that  the 
detection  statistics  are  the  output  of  the  flow  chart  in  fig,  Al, 


3l> 


Fig.  M(a-b),  Top-laval  flow  oharfc«  of  the  simulation 
programs*  CSOMCS  (a)  and  TABLE  (b). 
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START 


i 


Sntsr  fils  nsms 

Spaa if y  aampl 
number , of 
-and  range 

Lng  interval, 

samples 

covered 

Specify  type  of 
aperture  us ad 

ZI  Y 


Compute  ths  Tingle  integral 
ot  4w(t)  ualnfji^^ 


(iomputa  splint  cosmic  isn't  a" 
of  s  ( t)  and  i  ft) 

- - i  " 


Stors  samples  of  s©(t)  and 
positions  of  samples, 
ind  splins  oosfficisnts 
undsr  ths  given  file  nsms 


STOP 


Fig.  Ai(a-b),  Continued, 


41 


atois 

QMULT3 

QMUtTl 

IQHSCU 


it  oan  be  modified  easily  to  yield  the  estimation  statistics. 

A. 2  Quasi-Newton  Method 

There  are  several  versions  of  the  so  called  Quasi-Newton 
method  [8],  One  particular  version  has  been  written  by  Dr.  R.  w. 
Miller  in  the  Laboratory.  This  version  is  modified  here  to  meet 
the  requirements  of  the  present  simulation  study.  The  procedures 
are  as  follows! 

Step  It  Set  k«0i  read  in  the  initial  parameter  values ,  x°, 
and  the  number  of  targets ,  n. 

Step  2 1  Compute  J(x°)  and  VJ(x°)  according  to  Bqs.  (12), 
(15)  and  (18).  Find  H°  by  inverting  the  matrix 

i°(iij)  "  j|jrli“^xi' •  •  *  ,xi+^'  xi+l' ' "  ,x2n^’STx^‘^-0|]  ^'3"lr2,...,2n 

where  d  is  a  fixed  perturbation,  If  Q°  is  unin- 
vertable,  H°  is  chosen  as  a  diagonal  matrix  with 

H° (i,i)  -  d/| |?J(x°) | | 

where  |  j  *  |  |  is  the  Euclidean  norm. 

j# 

Step  3 i  Compute  the  increment  of  x  . 

&xk  -  Hk  J(xk) 


and  its- else, 


If  a  is  lass  than  a  preset  step  size  a ,  go  to 
step  4.  Otherwise,  multiply  Ax*  with  a/o  and  set 
o  ■  a. 

v 

Step  4 i  Compute  the  increment  of  J{x  ), 

AJ(xk)  -  Ax*  7J(x*). 

If  AJ(x*)  is  positive,  go  to  step  5.  Otherwise, 

"■  k 

replace  Ax  with 

v  7J(x*) 

AX*  -  0 - - - 

||7J(X*)|| 

where  $  is  a  preset  number. 

Step  5 i  Update  g, 

xk+1-x*+Ax* 
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Recompute  the  magnitude  of  the  increment r 


■  o  -  | | | | 

Step  7 1  compute  J(xk+1),  7J(xk,fl)  and 

k+i  .  (Aick+H\k)-  (Axk)T 
U*k)T  Hk^k 

provided  the  denominator  ia  not  equal  to  aero. 
Here  i 


»  7J(xk+1)  -  7J(xk). 

Step  8 i  if  e  <  cQ,  the  iteration  convergea  and  xk  ia  the 

deaired  eatimate.  Stop.  Otherwiae  proceed.  Note 
that  cQ  la  the  threahold  value  of  the  atop  criterion • 

step  9 1  if  k  <  k0,  aet  k  *  k+1  and  go  to-,  step  3.  other¬ 
wiae,  proaeed  to  the  following  step.  Here,  k0  ia 
a  preset  value  for  the  maximum  number  of  iterations. 

Step  10  s  Among  kQ  iterations  find  the  one  which  has  the 

greatest  value  of  J(jek),  Take  the  eatimate  of  this 
iteration  as  the  final  eatimate.  Then  atop. 

The  above  prooedurea  are  implemented  in  subroutines  QN1W  and 
ITBRAT.  The  neoeaaary  constants  are  preset  for  the  aimulation  as 
d  -  .003,  a  -  .5,  B  -  .13,  cQ  -  10“5,  and  k0  -  50. 
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A, 3  Computations  of  sQ(t)  and  *0(t) 

The  Quasi-Newton  procedure  requires,  at  eaoh  iteration,  the 
computations  of  J(x)  and  VJ(x)  which,  in  turn,  requires  those  of 
sQ (tg-T^)  and  s^tt^-T^)  for  fc-1,2, . *-,■  ,k  and *v*W  vIt^ir  •' 
very  time-consuming  to  directly  carry  out,  for  every  iteration, 
the  associated  single  and  double  integrals.  Although  the  computa¬ 
tion  time  can  be  reduced  significantly  by  using  the  high-speed 
convolution  method,  it  is  still  quite  noticeable.  An  alternative 
approach  is  to  compute  in  advanoe  samples  of  sQ(t)  and  sQ(t)  as 
wall  as  their  spline  coefficients,  and  store  them  in  files.  At 
the  beginning  of  the  simulation,  those  data  are  first  retrieved 
and  later  on,  whenever  needed,  and  «0 ( tA-r can  be 

computed  by  simple  interpolations  which  take  almost  no  time.  If 
the  number  of  stored  samples  is  large  enough,  the  interpolation 
would  provide  sufficient  accuracy, 

Since  s  (t)  and  e0(t)  would  be  only  computed  onoe  and  off-line 
the  computer  time  required  is  not  critical.  They  are  oomputed  by 
directly  carrying  out  the  double  and  single  intergrals  of  Bqs,  (3) 
and  (17}  using  the  Oausaian-Lagendre  quadrature  formula  [10,  11] , 
This  formula  gives  for  the  single  integral, 


X 


f! 


f (u)du 


(A.l) 
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the  following  approximation 


where  the  weights  A^n*  and  abscissas  x^n'  can  be  found  from  a 
standard  mathematical  table.  The  formula  is  exact  whenever  f (u) 
is  a  polynomial  of  degree  <  2n-l.  The  double  integral  of  the  form 

r  b  r  ♦<»») 

I  ■  /  /  f(u,v)dudv  (A. 3) 

J  b  (v) 


ia  approximated  by 


Here,  both  n  and  m  are  ohoeen  equal  to  16  for  the  current  applica¬ 
tion  »  which  is  shown  experimentally  to  be  sufficient.  Two  sub¬ 
routines  QMULT2  and  QMULT1  have  been  written  for  this  purpose. 

The  spline  coefficients  used  to  interpolate  the  set  of  pointa 
from  e0(t)  (or  a0(t))  are  computed  from  the  Quasi-Cubic  Hermits 
splines  [12 ].  The  cubic  spline  representing  the  function  between 
each  pair  of  given  points  is  determined  by  the  coordinates  and 
elopes  at  the  two  points.  The  slope  at  each  point  is  determined 


locally  by  the  point  in  gueition  and  two  point!  on  its  each  aide. 
The  resulting  curve  passes  through  all  the  given  points.  The 
subroutine  1QHSCU  which  is  available  in  the  IMSL  package  [13]  is 
employed  for  this  computation.  The  subroutine  SSDOT  called  by 
the  simulation  main  program  .performs  the  neoessary  interpolations 
to  yield  values  of  s^tt^)  (Eg.  (2))  and  3  s^t^/S  x^.  (Eq.  (16)), 
for  i»l,  2,..,,k  and  i«l,  2,,,.,2n  given  any  a  and  t.  Using  the 
output  of  SSDOT  the  subroutine  JFCN  computes  the  desired  J(x)  and 
VJ(«). 


A, 4  Choice  of  Initial  flusss 


The  subroutine  IC  is  employed  to  provide  the  initial  guess 
for  the  Quasi-Newton  procedure.  The  associated  principal  logic 
has  been  described  in  subsection  3,2,  In.  this  subsection,  the 
details  of  the  program  and  the  flow' chart  are  given. 

At  the  (i+1)  step  of  the  AIC  procedure,  the  initial  guesses 
of  the  first  i  targets  are  set  equal  to  the  estimates  from  the  itil 
step  and  the  initial  guess  for  the  (i+l)th  target  is  found  through 
the  pure  random  searoh.  From  Eq,  (12),  it  is  clear  that  max¬ 
imizing  J(x)  is  equivalent  to  minimising  the  following  function 


(|> 


i+1  - 

^Vo(VV 


(A.S) 


over  Sj  and  ,  j»l, , , , ,i+l. 
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APPENDIX  B i  PROGRAM  LISTINGS 

The  whole  eimulation  algorithm  ii  composed  of  two  main  pro¬ 
grams  *  TABLE  and  CSOMCS.  The  important  input  and  output  variables 
as  well  as  the  listings  of  theie  programs  are  given  in  this  appen¬ 
dix.  As  for  the  detailed  logic  of  the  algorithm*  readers  should 
refer  to  the  context  of  the  report.  The  comments  incorporated  in 
the  program  will  be  also  helpful  in  understanding  the  program. 

Important  Variables 


TABLE i 

IP 

■ 

indicator  of  the  type  of  aperture >  1  for  square 
aperture  and  2  for  annular  aperture 

E 

m 

obscuration  faotor  of  the  annular  aperture)  (D<E<1 

ANGLE 

m  . 

orientation  of  the  detector*  specified  by  the  angle 
(degree)  between  the  central. line  of  the  detector  and 
the  cross-scan  direction)  0<ANGLE<90. 

4SM  MS 

NPOINT 

m 

number  of  points  where  values  of  s0(t)  and  8Q(t)  are 
computed 

TDELTA 

m 

interval  between  a  pair  of  successive  points  in  units 
of  A/d 

TBEGXN 

m 

position  of  ths  first  point*  equal  to  (NPOINT-.I) 
•TDELTA/ 2 

ARG 

- 

array  containing  positions  of  the  NPOINT  points 

VALl 

m 

array  containing  the  corresponding  NPOINT  values 
of  s0(tt) 
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VAL2  ■  array  containing  the  NPOINT  values  of  s0(t^) 

SPLINI  =  array  of  dimension  NPOINTxS,  containing  the  spline 

coefficients  of  s^(t^)  .  v'  '••x  • 

SPLIN2  -  array  of  dimension  NPOXNTxl,  containing  the  spline 
coefficients,  of  s_  (t4 ) .  ; 

■  '0  X  r  •'  v.v  -v*  V’  ••  :  v  • 

......  .  ....  ; ...  ..  .  S- ':M:  r<-'  in;  ‘  1,-vJus  ty,* •  ■  v.As'JM.i.  v,«  •  V .  \.i  ;»u*. 1 

filen  ■  name  of  the  file  which  keeps  the  necessary  data 
(FILEN,  TSSGIN, ■ TDELTA,  NPOINT,  ARG,  VAL1, 
splinI  and  SPLIN2)  for  use  in  the  simulation  program 
CSOMCS 

CS0MC3 l 

FILEN,  TBEGIN,  TDELTA,  NPOINT,  ARG,  VAL1,  VAL2 ,  SPLINI, 

SPLIN2  as  defined  above  for  the  program  TABLE. 


The  following  parameter's  airs  specified.  ip  order  to  generate 
artificial  noisy  signal,.  .! :^y M  ' 

ix  »  seed  of  the  random  number  generator 

NLOOP  •»  number  of  Monte  Carlo  simulation  runs 
NTO  -  number  of  targets  present 

NS  -  number  of  samples 

DT  -  sampling  Interval 

AMP  «  array  containing  intensities  of  the  NTO  targets 


THETA  *  array  containing  positions  of  tha  NTO  targets 

# 

STD  "  standard  deviation  of  the  WON 

BT  »  position  of  the  first  sample,  equal  to  (NS-l)*DT/2 

'  The '  following  parameters  are  used  in  the  diteetioh/aitimatioh 
procedure i 

NT  -  number  of  targets  assumed 

N  -  total  number  of  free  parameters,  equal  to  2xNT 

.  \  .  •  ...  r.  *  . 

X  >  Vector  of  length  N,  The  intensity  and  position 

estimates  of  the  ith  target  are  kept  in  x(i]_  and. - 

X(i+NT).  - - - - 

NM _ _ _ _ «■ . preset  upper  bound  of  the  number  of  assumed  target 

.  •  '  i  .  ■ 

ETA  -  value  of  the  aic  penalty  constant 

AIC  ■  Value  of  the  Akaike  information  Criterion 

SSQ  *  value  of  J(x)  in  Bq,  (12) 

SSQN  »  value  of  SSQ  under  the  assumption  that  no  target  pre- 

■  sent 

S  *  vector  of  length  NS  containing  values  of  s^(t)  taken 

in  the  N8  points 
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SDOT 


XBAR 

XVAR 

1COUNT 

EPS 

LIMIT 

ITRL 

Listings 


vector  of  length  NSxN  containing  partial  derivatives 
of  sd (t) ; 

3sd(tl)  3aa(tt)  3sd(tt> 


3W 


3a. 


9  •  •  •  9 


>aNT 


3t. 


t  •  •  •  t 


are 


at 


NT 


stored  sequentially  in  the  locations  beginning  at 
SDOT(  (A-l)  *N+1) 

=  vector  containing  mean  values  of  the  estimates  in  all 
the  models  assumed 

=  vector  containing  variances  of  the  estimates 

i.L 

=  array  whose  i  element  containing  the  probability  of 
identifying  i  targets 

=  smallest  increment  used  in  convergence  test  (for 
QNEW  subroutine) 

«  allowed  largest  number  of  iterations  (for  QNEW  sub¬ 
routine) 


number  of  the  random  searches  used  in  subroution 
IC  for  choosing  the  initial  guess 
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MAIN  TABLE 


C 

C* 

c 

C  ...  DESCRIPTION! 

C  THIS  PROGRAM  COMPUTES  AND  8T0RE8  THE  OPTICAL  PUL8E  8HAPE  AND 
C  ITS  FIRST  DERIVATIVE  U.R.T.  ANGULAR  POSITION  FOR  FUTURE  U8E. 

C  TilE  SINGLE  INTEGRAL  OF  EG. (17)  IS  EVALUATED  BY  USING  QAUB8- 
C  LEGENDRE  FORMULA.  SUBROUTINE  GL016  PROVIDES  A  TABLE  OF  THE  16- 
C  POINT  GAUSS-LEGENDRE  FORMULA. 

C  THE  DOUBLE  INTEGRAL  OF  EG. (3)  18  COMPUTED  USING  THE  SAME  METHOD. 

C 

EXTERNAL  FCN.FUP.FLO.FCT 
DOUBLE  PRECISION  DX< 16) .DA< 16) 

DIMENSION  X  <16  >  >  A<16).  FILEN<3) 

DIMENSION  VAL1 (1024 ).VAL2 <1024 >»ARGd024)»SPLINl <1024. 3) »8PLIN2< 
11024.3) 

COMMON  /PSF/  IP.E  /DETOR/  CENTER. ANGLE. BETAX 
URITE(6. 10) 

10  F0RMAT(//10X. '***  TABLE  ***'//> 

URITE<6.20> 

20  FORMAT ( 1 X . ' READ  IN  FILE  NAME') 

READ(5.30)  FILEN 
30  FORMAT (3A4) 

URITE(6»40)  FILEN 
40  FORMAT ( IX. 'FILE  NAME-'.3A4> 

TDELTA  -  .025 
NP0INT-1024 

RANGE  »  <NP0INT-1 )  *  TDELTA 
T BEGIN  -  -RANGE  /2. 

BETAX  -  2. 

BETAY  -  6. 

BETAYH  -  BETAY/2. 

BETAYL  -  -BETAYH 
URITE<6»50) 

50  FORMAT! IX. 'READ  IN  IP!  1  FOR  RECTANGULAR.  2  FOR  ANNULAR  ') 

READ  (5.4i)  IP 
URITE  <6.60 )  IP 
60  FORMAT  (IX, 'IP-  '.12) 

IF  (IP  .EO.  1  )  GO  TO  80 
URITE (6. 65) 

65  FORMAT! IX. '  READ  IN  0B8CURATI0N  FACTOR  ') 

READ<5.*>  E 
URITE(6»70>  E 

70  FORMAT (IX. '  E-  '»F4.2) 

80  URITE(6.90) 

90  FORMAT (IX. 'READ  IN  DETECTOR  ANGLE< DEGREE) ' ) 

READ (5.*)  ANGLE 
URITE(6» 100)  ANGLE 
100  F0RMAT<1X» 'ANGLE"' .F6.2) 

C 

C  OBTAIN  THE  GAU88-LEGENDRE  UEIGHT8 
C 

CALL  GL016(DX.DA.-1 .DO.  l.DO  ) 

DO  110  1-1.16 
X(I>  -  DX< I ) 

110  A(I)  •  DA(I) 

MM  -  16 
C 

C  COMPUTE  THE  PUL8E  8HAPE.VAL1  AND  ITS  DERIVATIVE  VAL2 
C 

DO  120  I-1.NP01NT 

CENTER  -  TBEOIN  +  TDELTA*  <I-1) 
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ARG < I )  «  CENTER 

MAI.  1  <  1 )  -  0MULT2 <  FCN , BETA YL , BETA YH » PUP .FLO.X.A.MM) 

MAL2<  l >  -  QMULT1  (FCT ,  BETAYL.BETAYH.X, A.MM) 

120  CONTINUE 
C 

C  COMPUTE  THE  SrLINE  COEFFICIENT  OF  MALI.  SUB  IQHSCU  IS  IN  IMSL. 

C 

CALL  IOHSCUIAROf MALI .NPOINT .3PLIN1 .NPOINT » IER) 

C 

C  COMPUTE  THE  SPLINE  COEFFICIENT  OF  MAL2. 

C 

CALL  I GHSCU < ARG »MAL2» NPOINT .SPLIN2. NPOINT » IER) 

C 

C  WRITE  OUT  THE  DATA 
C 

WRITE (6. 130)  TBEOINf TDELT A* NPOINT 
130  F0RMAT</lX,2E14.7.I8/> 

DO  150  1=1 .NPOINT 

WRITE (6. 140)  I ?ARG< I ) .MALI <1 >  »MAL2( I > » (SPLIN1 ( I* J) . J=1 .3) . (SPLIN21 
*I,J) , J= 1 . 3 ) 

140  FORMAT (1X,I5,2X.9E12<5) 

150  CONTINUE 
C 

C  STORE  DATA  UNDER  THE  GIMEN  FILE  NAME 
C 

WRITE  <  8. 1 AO )  FILEN.  TBEG IN .TDELTA. NPOINT 
DO  170  1=1 .NPOINT 

WRITE (8. 180)  ARG< I ) .MALI < I ) .MAL2< I  >  »  < SPLIN1 ( I » J) ,J=1.3> . 

1  <  SPLIN2  <  I »  J) , J=1 , 3 ) 

170  CONTINUE 

160  FORMAT (3A4.2E14.7.I6) 

180  FORMAT  <9E14.7) 

STOP 

END 


C*  FUNCTION  QMULT2 

C 

C  COMPUTE  THE  DOUBLE  INTEGRAL  SS  F(X.Y>  DXDY.  AA.LE. Y.LE.BB, 
C  FL<Y) .LE.X.LE.FU(Y) 

C 

FUNCTION  QMULT2 < FCN . AA , BB . FU . FL . X . A . MM > 

DIMENSION  X(1).A(1) 

HI  =  <BB-AA)/2. 

G1  «  ( BB+AA ) /2 . 

01  =  0, 

DO  4  1=1, MM 

UI  =  H1*X<I)  +  G1 

AI  »  H1*A( I ) 

D  »  FU<UI ) 

C  «  FL<UI) 

H  =  <D-C>/2. 

G  «  (D+O/2, 

0  =  0. 

DO  2  J-l.MM 
MJ  «  H*X<J)  +  G 
2  0  =  0  +  A < J )  *  FCN(MJ.UI) 

4  01  «  01  +  AI*H*0 

0MULT2  »  01 
RETURN 
END 
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o  no  o  non 


FUNCTION  FCN 


C 

c 


10 


32 

36 


4  2 
A6 

20 


50 


COMPUTE  F<X» Y)  WHICH  IS  USED  IN  QHULT2 

FUNCTION  FCN<XfY) 

DOUBLE  PRECISION  MMBSJ1 r DR 
COMMON  /PSF/  IPfE 
SCALE  »  3.141593 
00  TO  < 10*20) » IP 
IF  <X  .EQ.  0.)  GO  TO  32 
XI  =  X  *  SCALE 
FX  =  <SIN<X1)/X1)**2 
GO  TO  36 
FX=1 . 

IF  (Y  .EQ.  0. )  GO  TO  42 
Y1  =  Y*  SCALE 
FY  »  (SIN<Y1  >/Yl>**2 
GO  TO  46 
FY  =1. 

FCN  -  FX*FY 
RETURN 

R  =  SORT <  X**2  +  Y**2  )  *SCALE 
IF  < R  .NE.  0.)  GO  TO  50 
FCN*=1 . 

RETURN 
DR-R 


MMBSJ1  COMPUTES  THE  BESSEL  FUNCTION  OF  FIRST  KIND r EXISTING  IN  IMSL. 

BJ1=MMBSJ1 (DR» IER1 ) 

IF< E. EQ. 0 . )  GO  TO  60 
DR=R*E 

B J2=MMBSJ1 <  DR» IER2 ) 

GO  TO  70 
60  IER2-0 

BJ2  ■  0. 

70  IF  (IER1  .NE.  0  .OR.  IER2  . NE .  0)  GO  TO  BO 
TEMP  ■  2,/((l.-E**2)*R) 

FCN  -  < ( BJ1  -  E*BJ2>  *TEHP>**2 
RETURN 

BO  WRITE  <  6»  90) 

90  FORMAT  < //10X  t '  SUB  MMBSJ1  ERROR  '> 

CALL  EXIT 
END 


*  FUNCTION  FUP 

THE  UPPER  BOUND  USED  IN  QMULT2 
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FUNCTION  FUP<Y) 

COMMON  /DETOR/  CENTER  t ANGLE  *  BETAX 
FUP  -  CENTER  +  BETAX/2 . 

IF < ANGLE. EQ.O. )  GO  TO  10 
FUP-FUP+ Y*TAN ( ANGLE* .01745329  > 
RETURN 
END 
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C*  FUNCTION  FLO 

C 

C  THE  LOWER  BOUND  USED  IN  QMULT2 
C 

FUNCTION  FLO(Y) 

COMMON  /DETOR/  CENTER  t ANGLE t BET AX 
FLO  =  CENTER  -  BETAX/2. 

IF < ANGLE. EO. 0. )  GO  TO  20 
FLO=FLO+Y*TAN< ANGLE*. 01 745329) 

20  RETURN 
END 


C*  FUNCTION  QMULT1 

C 

C  COMPUTE  THE  SINGLE  INTEGRAL  S  FCY)  DY 
C 

FUNCTION  GMULT1 < FCT f AAf BB f X  f A f MM) 
DIMENSION  X<1)fA<1> 

Hl=<BB-AA)/2. 

Gl=<BB+AA)/2. 

01=0. 

DO  4  1=1 »MM 

UI=H1*X<I)+G1 

AI=H1*A<I> 

4  01=01+AI*FCT<UI) 

QMULT1=Q1 

RETURN 

END 


C*  FUNCTION  FCT 

C 

C  COMPUTE  F<Y)  WHICH  IS  REQUIRED  BY  QMULT1 
C 

FUNCTION  FCT  <  Y  ) 

COMMON  /DETOR/  CENTER  t ANGLE f  BETAX 
X=CENTER 

IF(ANGLE.E0.0. )  GO  TO  30 
X  =  X  +  Y*TAN<ANGLE  *  .0174532?) 

30  XI  =  X  +  BETAX/2. 

X2  =  X  -  BETAX/2. 

FCT  =  FCN(XIfY)  -  FCN<X2fY> 

RETURN 

END 


C*  SUBROUTINE  GL016 

C 

C  PREPARE  COEFFICIENTS  OF  THE  16-POINT  GAUSS-LEGENDRE  FORMULA 
C 

SUBROUTINE  GL016 < X f A f Cf D) 

DOUBLE  PRECISION  C . D f X < 1 > » A< 1 > f XX ( 8 ) f AA< 8 > 

DATA  XX/ 

t  .98940093499 1649932596 1541 734D0  f 
*  • 9445750230732325760779884 1S5D0  f 
t  • 863631 20238783 17438 8046 78977D0  t 
%  . 755404408355003033895101 194 8D0  f 
i  . 61 7876244402643748446671 7640D0  f 
%  . 45801 677765722738634241 944 29D0  f 
t  . 20 16035507792589 1 323046050 14D0  f 
t  . 9501 25098376374401 853 1933542D-1  / 
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DATA  AA/ 

t  .271 524594 11 754094851 78037245D-1* 
l  • 62253 5 2393 864 78928628 4 383699D-1 t 
t  •9515851 1682492784B0992510760D-1 t 
%  . 124A289712555338720524762821D0  t 
S  .  14959598881<S5747320813017305D0  » 
t  . 1691565193950025381893120790DO  * 

S  . 1826034150449235888667636479D0  t 
S  • 18945061045506B4962853967232D0  / 

DMC  -  • 5D0#  <  D-C ) 

DPC  =>  •  5D0*  ( D+C ) 

DO  2  1=1 t 8 
NI  =  17-1 

X<I)  =  -DMC*XX<I)  +  DPC 
X  <  NI )  =  DMC#XX < I )  +  DPC 
A<I)  =  DMC*AA<I) 

A  <NI >  =  DMC*AA<I> 

RETURN 

END 
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MAIN  CSOHCS 


C 

C* 

C 

C  -  DESCRIPTION  t 

C  THIS  PROORAM  IS  A  MONTE-CARLO  SIMULATION  OF  THE  INTENSITY  AND  ANGU 
C  LAR  LOCATION  ESTIMATION  OF  CLOSELY  SPACED  OPTICAL  TARGET8. 

C  A  SET  OF  ARTIFICIAL  OBSERVATIONS  ARE  FIRST  GENERATED  FOR  A  GIVEN 

C  CONFIGURATION  (NO.  OF  TARGETS.  INTENSITIES  AND  LOCATIONS). 

C  USING  THESE  DATA  AND  ASSUMING  NO.  OF  CSO'S.  A  GAUSI-NEWTON  ALGORI 
C  THM  IS  EMPLOYED  TO  SEARCH  FOR  THE  PARAMETER  SET  WHICH  MAXIMUMS  THE 
C  LIKELIHOOD  FUNCTION. 

C  THE  MEANS  AND  VARIANCES  OF  THE  ESTIMATES  ARE  ALSO  COMPUTED. 

C  AKAIKE  INFORMATION  CRITERION  IS  COMPUTED  FOR  VARIOUS  MODELS. 

C 

C 

INTEGER  1COUNT  <  5 ) 

REAL*4  AMP <4  > » THETA <4 ) . X ( 8 > » S<64 > . SDOT (512) 

REAL*4  XBAR< 20 )  .XVAR<20>  . AML (100 . 5 )  . DUM< B )  » THUO )  »FILEN<3) 

COMMON  /DATA/Y<64) .NS.STD  /SAMPLE/BT.DT  /WORK/SDOT 
COMMON  /SPLINE/  TBEG IN » TDELTA . NPOINT . ARG ( 1024 > » VAL 1 ( 1024 ) » VAL2 ( 102 
1  4 ) .SPLIN1 ( 1024.3) .SPLIN2< 1024. 3  > 

WRITE<6» 10) 

WRITE  (6.20) 

10  FORMAT </5X. '*****  RESULTS  FROM  CSOMCS  #****') 

20  FORMAT < 10X. 'CSO  MONTE  CARLO  SIMULATION') 

WRITE<6.30> 

30  FORMAT ( 14X. 'WGN  <STD  UNKNOWN)') 

C 

C  READ  IN  STANDARD  PULSE  SHAPE  FROM  THE  FILE  WHICH  HAS  BEEN  CREATED  BY 
C  PROGRAM  TABLE. 

C 

READ  <8.50)  F I LEN , TBEG I N . TDELTA . NPO I NT 
DO  40  1=1. NPOINT 

READ (8. 60)  ARG( I ) , VAL1 ( I ) . VAL2 ( I ) »  < SPLIN1 < I »  J) .J-1.3) . 

1  <  SPLIN2C IfJ)»J=1.3> 

40  CONTINUE 

50  FORMAT (3A4.2E14.7. 16) 

60  FORMAT ( 9E 1 4 . 7 ) 

WRITE<6»70)  FILEN 

70  F0RMAT</12X. 'FILE  NAME) ' »3A4/) 

C  WRITE<6»*)  TBEGIN. TDELTA. NPOINT 

C  DO  12  I»l. NPOINT 

C  WRITE<6. 13)  T  » ARG< I ) . VAL1 ( I > . VAL2 ( I >  » (SPLIN1 <I.J)»J*1»3). 

C  *<SPLIN2<I.J). J-1.3) 

C13  F0RMAT<1X»I5.9E12.4) 

C12  CONTINUE 

WRITE(t..80> 

80  FORMAT < lX'ENTER  THE  PENALTY  CONSTANT  OF  AIC') 

READ<5»*>  ETA 
WRITE(6»90)  ETA 
90  FORMAT (2X.F5. 2) 

WRITE(6. 140) 

140  FORMAT ('  ENTER  SEED  FOR  RANDOM  NUMBER  GENERATOR') 

READ(SrB)  IX 
WRITE<6.1S0)  IX 
150  FORMAT  <3X» 19) 

WRITEC6.160) 

160  FORMAT < IX. 'ENTER  NO.  OF  MONTE-CARLO  LOOPB') 

READ(5.*)  NLOOP 
WRITE<6. 170)  NLOOP 
170  FORMAT <3X» 14) 
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WRITE!6r 180) 

180  FORMAT < '  ENTER  DATA  TO  GENERATE  ARTIFICIAL  SIGNAL! '/ 

1  5X* 'NT »N8»DELTAr AMP!NT> »THETA!NT>  »8TD' > 

READ!5»*>  NT0»N8»DT*  <AMP< I ) » I-l # NTO) » < THETA < I ) » I-l »NTO) *8TD 
WRITE!6»190)  NTO»NS»DT » ! AMP! I ) » 1*1 »NTO> »  ! THETA! I ) » l-l » NTO) » STD 
190  FORMAT ! 1X» 13. I5»9F8.2> 

BT— !NS-l)*DT/2. 

NM-NTO+1 

NH2-0 

DO  200  I-1»NM 
200  NM2-NH2+IB2 

DO  210  I-1.NM2 
XBAR!I>-0. 

210  XVAR!I>-0. 

NM1-NM+1 
DO  220  I-lfNHl 
220  ICOUNT ( I )-0 
EPS-1. E-5 
LIMIT-50 
ITRL-50 
C 

C  SIMULATE  NOISELESS  OBSERVATIONS 
C 

CALL  SSDOT !NTO» AMP . THET A » 8 » SDOT » N8 » 0 ) 

WRITE!6>230) 

230  FORMAT!/'  ESTIMATION  STARTS!  ') 

DO  400  NL-lrNLOOP 
WRITE ! 6»  240  >  NL.IX 

240  FORMAT  < / 1 X  t ' NLOOP- ' . 1 4 » 5X  f ' IX- ' . 1 10  > 

C 

C  SIMULATE  ARTIFICIAL  NOISY  DATA 
C 

DO  250  I-l. NS 
SMEAN  -  8(1) 

CALL  GAUSS!  IX . STD . SMEAN » Y ! I ) ) 

250  CONTINUE 
C 

C  COMPUTE  880  AND  AIC  FOR  NT-0 
C 

S8QN-0. 

DO  260  I-l. NS 
260  SSQN-88QN-Y ! I ) **2 

C 

C  VARIANCE  OF  THE  NOISE  IS  UNKNOWN 
C 

A I C-NB$ ALOG ! -SSQN/N8 ) 

C 

C  "FIANCE  OF  THE  N0I8E  IS  KNOWN 
C  AIC— 8S0N/3TDM2 

C 

WRITE!6»270>  SSQN.AIC 

270  F0RMAT!1X»'880N-'.E15,6»3X» 'AIC-'»E15.6> 

C 

C  APPLY  THE  AIC  PROCEDURE 
C 

NT-1 

280  N-NT42 

WRITE!6f 290)  NT 

290  FORMAT! 12 .'-TARGET  MODEL!') 

CALL  IC!X.N.88Q. IX. ITRL) 

WRITE! 6. 300)  (  i).I-l»N> 

300  F0RMAT!1X»'IC!  .8E12.4) 
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c 

C  USE  BM2  FOR  A  SO  THAT  B  18  NOT  CONSTRAINED  TO  BE  P08XTIVE 

C  DO  310  I-1*NT 
310  X<I)-SQRT<X<I>> 

CALL  QNEUT  <X»N*S80»EPS*LINIT. ITER.IER) 

C 

C  RESTORE  A-BM2 
C 

DO  320  I-l.NT 
320  X(I)-X(1)M2 

CALL  ORDER(X*N> 

WRITE (6 *330)  <X(I)*I-1*N) 

330  FORNAT <2X* 'ESTIMATED  PARAMETERS! ' *8E12.4> 

C 

C  VARIANCE  IS  UNKNOWN 
C 

AICNT  -NSt ALOG  (  -S8Q/N8  )  +  ETA<I<2<NT ) 

C 

C  VARIANCE  IS  KNOWN 
C  AICNT— SS0/8TD«2  +  ETA<<2<NT> 

C 

WRITE (6*340)  SSQ* AICNT *  ITER 

340  FORMAT (  2X  *  ' S8Q-  ' *E14.6*5X* 'AIC-' *E14.6*5X* ' ITER.-' » 14) 
IF (AICNT  .GE.  AIC)  GOTO  360 
AIC-AICNT 
DO  350  I-1*N 
350  DUM<I)-X(I) 

NT-NT+1 

IF  (NT  .GT.  NM)  GO  TO  360 
GO  TO  260 
C 

C  UPDATE  THE  COUNTER 
C 

360  I COUNT (NT)  -  ICOUNT(NT)  +1 

IF (NT  .EQ.  1)  GO  TO  400 
C 

C  KEEP  DATA  FOR  COMPUTING  8AMPLE  MEAN  AND  VARIANCE  OF  THE 
C  ESTIMATED  PARAMETERS 
C 

K-0 

DO  370  I«2*NT 
370  K-K-KI-2X2 
N»(NT-1><2 
DO  380  I**l  *N 

XBAR (K+I >*XBAR(K+I )+DUM< I ) 

380  XVAR (K+I ) -XVAR (K+I  )+DUM(  I  )«2 
400  CONTINUE 
C 

C  COMPUTE  AND  OUTPUT  SIMULATION  8TATX8TXC8 
C 

WRITE(6*410) 

410  F0RMAT(///1X'8TATI8TIC8' ) 

NT-0 

WRITE(6X20> 

420  FORMAT <2X* 'O-TARGRT  MODEL!') 

PCT- 100 . BICOUNT < 1 ) /FLOAT  <  NLOOP ) 

WRITE(6*440)  NT*PCT 
K-0 

DO  490  NT-1*NM 
N-NTB2 

WRITE(6*430>  NT 

430  F0RMAT(/1X*I2* '-TARGET  MODEL!') 
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COUNT-FLOAT ( IC0UNT(NT+1  >  ) 

PCT-COUNT/NLOOP* \ 00 . 

WRITE(6t 440)  NT »PCT 

440  FORMAT <3X» 'PR0B( ' » II t ' )-' »F6.2» > 

IF<PCT .EQ.O. )  00  TO  480 

DO  450  1-1 »N 

Il-K+I 

XBAR (11) -XBAR (ID /COUNT 

XVAR ( 1 1 ) - ( XVAR (ID -COUNTBXBAR  < Il)**2) /COUNT 
450  XUAR(  II  )-8QRT (XVAR(  ID  ) 

WRITE(6.460)  (XBAR(K+I >»I«1*N> 

460  FORMAT (3X» 'MEAN  i ' »8( 1X>E12.4) > 

WRITE ( 6  *  470 )  ( XVAR ( K+I ) » 1-1 »N ) 

470  FORMAT ( 3X» 'STD $  ' 1 0( 1X*E12.4) ) 

480  K-K+N 
490  CONTINUE 
CALL  EXIT 
END 


C*  SUBROUTINE  QNEU 

C 

SUBROUTINE  QNEUT ( X » N » X J » EPS » L I M I T » I CONT » I ER > 

THIS  SUBROUTINE  IMPLEMENTS  QUA8I-NEWT0N  METHOD  TO  FIND  A  MAXIMUM  OF 
A  FUNCTION. 

Xt  INITIAL  QUE8P( INPUT) ♦  LOCATION  OF  MAX I MUM (OUTPUT) 

Nt  DIMENSION  OF  X 
XJJ  VALUE  OF  MAXIMUM  'OUTPUT) 

EPS J  SMALLEST  INCREMENT  USED  IN  CONVERGENCE  TEST 
LIMIT!  MAXIMUM  NUMBER  OF  ITERATIONS 
ICONT!  NUMBER  OF  ITERATIONS 
IER!  ERROR  CODE 

REQUIRES  A  SUBROUTINE  JFCN(X'NfXJtDJ)  WHERE 
X!  ANY  POINT  IN  THE  PARAMETER  8PACE 
N!  DIMENSION  OF  X  AND  DJ 

xj:  value  of  the  function  at  x  (output) 

DJ!  GRADIENT  EVALUATED  AT  X  (OUTPUT) 

DIMENSION  X(l)fDJ(8) »DJ1 (8) »X1 (8) 

DIMENSION  B<8»8>rBI(8F8>fFJ(51)fFX(51r8) 
C0MHQN/NAIN1/NDIM/IN0UT/KIN»K0UT 
COMMON/PARAM/S  »D8*DSHfI0UT 

MAX(N>-8»  FJ(LIMIT+1 >  »  FX(LIMIT+1 »N> 

/MAIN1/  AND  /INOUT/  ARE  LINKED  TO  GMINV  SUB. 

NDIM  -  8 
KIN  -  5 
KOUT  -  6 
C 

C  FOLLOWING  C0N8TANT8  ARE  U8ED  TO  CONTROL  NUMERICAL  PROCEDURES  OF  QNEW 
C  S!  INITIALIZATION  PERTURBATION 

C  DSt  INCREMENT  U8ED  WHEN  NOT  NEAR  A  MAXIMUM 

C  DBM!  LARGE8T  INCREMENT  ALLOWED 

C  IOUT-PRINTOUT  CONTROL 

C 

I OUT-1 
8- <005 
D8-.15 
D8M-.5 
IER  -  0 
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DO  10  I-l.N 
FX(1.I)-X(I> 

CALL  JFCN(X.N.XJO.DJ) 

FJ(1)-XJ0 

IF( IOUT . GE . 1 >  URITE(6.20>  XJO 
20  FORMAT ( IX. 'SSQ0-' .E13>6) 

C 

C  INITIALIZE  D 
C 

DO  40  I-l.N 
DO  30  J-l.N 
30  X1<J>  -  X<J) 

X1<I>  -  X<I>  +  S 
CALL  JFCN(X1»N»XJ»DJ1> 

DO  40  J-l »N 

40  BI( J. I)  ■  (DJI ( J>-DJ<J) )/8 

CALL  GMINV(N.N.BI.B.MR.l) 
if<mr.lt.n;  go  to  60 
DO  50  I-l.N 
DO  50  J-l.N 
50  B(I»J)  »  -B(IfJ) 

GO  TO  100 
C 

C  ALTERNATE  INITIALIZATION 
C 

60  DJM  -  0. 

DO  70  I-l.N 

70  DJM  -  DJM  +  DJ(I>**2 

DJM  -  8/SQRT (DJM) 

DO  90  I-l.N 
DO  80  J-lrN 
80  B(I.J)  -  0. 

90  B(I.I)  -  DJM 

C 

C  ITERATE  TO  SOLUTION 
C 

100  DO  130  ICONT-1. LIMIT 

CALL  ITERAT(X.B.DJ. XJr DXN.N.IR) 

IC0NT1-IC0NT+1 
FJ( ICQNT1 )-XJ 
DO  110  I-l.N 
110  FX( IC0NT1 .1 )-X( I ) 

IF(I0UT.0E.3)  URITE(6» 120)  ICOUNT.XJ. (X( I > .1-1 .N) 
120  FORMAT ( 1X.I4. 2E12«4.4E12«4/(29X.4E12»4) ) 
IF(IR.NE.O)  80  TO  150 
IF (DXN.LT »EP8>  GO  TO  190 
130  CONTINUE 

IF(IOUT.GE.l)  WRITE (6. 140)  LIMIT 
IER  -  2 

140  FORMAT (2X. 'NO  CONVERGENCE  IN'. 14.'  ITERATIONS') 

GO  TO  160 
150  IER-3 

160  XJ-FJ(l) 

KMAX-1 

DO  170  1-2. IC0NT1 
IF(FJd).LE.XJ)  GO  TO  170 
XJ-FJ(I) 

KMAX-I 

170  CONTINUE 

DO  180  I-l.N 
180  X ( I ) -FX ( KMAX . I ) 

RETURN 
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190  IF(XJ.GE.XJO)  RETURN 
IER  -  1 

IF ( I OUT • GE • 1 >  WRITE (6f 200) 

200  FORMAT < IX 'NOT  GLOBAL  MAXIMUM') 

RETURN 
END 


C*  SUBROUTINE  ITERAT 

C 

C  CALLED  BY  QNEU 
C 

SUBROUTINE  ITERAT<Xf Bf FfXJ.DXNfNfIR) 

DIMENSION  X(l>*B(8»8>»F(l)fDF(8>»DX(8)>BF(8>»DXB(8) 

COMMON/PARAM/S . DS  f  D8M » IOUT  /8AMPLE/BT 

IR-0 

DXN  ■  0. 

DO  20  I-lrN 
DX< I )  -  0. 

DO  10  J-1»N 

10  DX ( I >  -  DX< I )  +  B(IfJ)*F(J> 

20  DXN  -  DXN  +  DX(I)**2 
DXN  -  SORT (DXN) 

IF(DXN.LE.DSM)  00  TO  SO 
C 

C  INCREMENT  TOO  BIG 
C 

DO  30  I-1»N 

30  DX< I )  -  DX< I ) tDSM/DXN 

DXN  -  DSM 

IF (TOUT ,0E.3>  WRITE(At 40) 

40  FORMAT < '  INCREMENT  TOO  BIO') 

C 

C  CHECK  FOR  CORRECT  MOTION 
C 

50  XDF  -  0. 

DO  AO  I-l.N 

AO  XDF  =  XDF  +  DX(I)*F(I> 

IF (XDF .GT .0. )  00  TO  100 
C 

C  MOTION  IN  WRONG  DIRECTION)  CHANGE  TO  THE  DIRECTION  OF  GRADIENT 

C 

IF ( IOUT  »GE»3>  WRITE( Af 70) 

70  FORMAT ( '  NEAR  MINIMUM') 

FM  -  0, 

DO  80  I-1'N 

80  FM  -  FM  +  F ( 1)4*2 

DXN  -  DS 
FM  -  D8/8QRT(FM> 

DO  90  I-1»N 
90  DX( I )  -  FM*F(I> 

100  DO  110  1*1 *N 

110  X(I)  -  X(I>  +  DX( I > 

C 

C  CONSTRAINT  THE  TARGET  P08ITI0NS  WITHIN  THE  RANGE  (BTf-BT) 

C 

NT1-N/2+1 
I FLAG-0 

DO  120  I-NTlrN 
AB8X-AB8(X(I)) 
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XF<ABSX.LE.-BT>  00  TO  120 
IFLAG-1 

SIGNX“X< I )/AB8X 

TEMP— SIGNXBBT 

DX< I )«DX< I )+TEMP-X< I ) 

X<I)-TEMP 
120  CONTINUE 

IF<IFLAG.EQ.O>  00  TO  ISO 
DXN-O. 

DO  130  I-lfN 
130  DXN-DXN+DX<I>**2 

IF<DXN.NE.O. )  00  TO  140 
IR-1 

XJ—  1.E75 
RETURN 

140  DXN«SGRT<DXN> 

150  CALL  JFCN  <  X  f  N  f  X J  f  DF  > 

DO  160  1*1 fN 

DF< I )  -  DF( I )  -  F< I > 

160  F( I )  -  DF( I >  +  F< I ) 

DO  170  1-1  fN 
BF< I )  -  0. 

DXB<  I )  «*  0. 

DO  170  J*1  fN 

BF < I )  -  BF ( I )  +  B(I»J>*DF(J) 

170  DXB< I >  -  DXB(I)  +  DX<J)*B<J'X) 

A1  -  0. 

DO  180  I»1>N 

180  A1  -  A1  +  DX<I)*BF<I> 

IF< A1 .NE iOi)  00  TO  190 

IR-2 

RETURN 

190  A1  -  1./A1 

DO  200  1*1 fN 
DO  200  J-1fN 

200  B(IfJ)  -  B< I f J>  -  A1*(BF( I )+DX( I ) >*DXB( J) 

IF < IOUT . GE. 4 )  URITE<6f210)  ( <B( Jf X ) fX-1 fN) f J-l fN) 
210  FORMAT ( 2E20 • 4 ) 

RETURN 
END 


C*  SUBROUTINE  GAUSS 

C  GENERATE  GAUSSIAN  DISJE^uflON  RANDOM  NOISE 
SUBROUTiHS-USuSSf IXfSf AMfM) 

10  1*1  f  12 

CALL  RANDU<IXfIYfY> 

IX-IY 

10  A  -  A  +  Y 

V  -  <A-6.0>  *8  +AM 

RETURN 

END 
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c*  SUBROUTINE  JFCN 

C 

C  COMPUTE  J<X>  AND  ITS  GRADIENTS 

C 

C 

SUBROUTINE  JPCN  <X»N»880»DJ> 

DIMENSION  AMP<4) » THETA<  4 ) » X ( 1 ) » DJ ( 1 ) *  S ( 64 ) r BOOT (312) 

COMMON  /DATA/  Y<64)»NS»STD  /W0RK/8D0T»8 

NT-N/2 

DO  10  K-1»NT 
AMP<K)«X<K>**2 
10  THETA<K>-X<K+NT> 

CALL  SSDOT  <NT , AMP t THETA f St SDOT.NSf 1 > 

SSO-O. 

DO  20  JL-l.NS 
DHL  »  Y < JL )  -S< JL> 

SSO-  SSG-BML**2 
20  CONTINUE 

DO  40  IR=1 t N 
SUM-0. 

DO  30  JL-1 »NS 
NN=N*(JL-1>  +  IR 
BML  -  Y < JL )  -S(JL) 

SUM-SUM+SDOT ( NN ) *BML 
30  CONTINUE 

C 

C  VARIANCE  IS  UNKNOWN 
C 

DJ<  IR)«*-SUM/SSQ#NS 
C 

C  VARIANCE  IS  KNOWN 
C 

C  DJ(IR) 

40  CONTINUE 

RPTURi r 
' END 


C*  SUBROUTINE  ORDER 

C 

C  ARRANGE  THE  ORDER  OF  THE  TARGETS  ACCORDING  TO  THEIR  IN-SCAN  POSITIONS 
C 

SUBROUTINE  ORDER<X*N> 

DIMENSION  X<1) 

NT-N/2 

IF(NT.LE.l)  RETURN 

I 1-NT+l 

I2=»N-1 

DO  20  I-Il » 12 
Jl-I+1 

DO  10  J-J1,N 

IF<X<  J) »LE.X< I )  )  GO  TO  10 
TEMP-X( I > 

X<I)=X<J> 

X< J)-TEMP 
TEMP-X(I-NT) 

X(I-NT>-X< J-NT) 

X< J-NTJ-TEMP 
10  CONTINUE 

20  CONTINUE 

RETURN 
END 
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SUBROUTINE  IC 


C* 

H  PROVIDE  INITIAL  OUESS  FOR  QNEW 
C 

SUBROUTINE  IC<X,N,8SQ, IX, ITRL) 

DIMENSION  X<1),S<64>,SD0T<512>,Y1<64> 

COMMON  /SAMPLE/BT  /W0RK/8D0T»S  /DATA/Y(64) ,  NS 

NT-N/2 

AL-BT 

AU— BT 

IF<NT.NE.1>  00  TO  20 
XJO-1 .E75 
DO  10  I=1*NS 
10  Y1 ( I )=Y< I ) 

GO  TO  50 
20  XJO=-SSO 

NT1=NT-1 

CALL  SSDOT (NT1,X<1),X< NT ) iS> SDOT , NS , 0 ) 

DO  30  1=1 » NS 
30  Y1 ( I )=Y< I )-S( I ) 

DO  40  1  =  1  »NT1  _ _ 

40 _ _JLLMr-£-f-*rX ) 

IFLAG=0 

DO  90  K=1 » ITRL 
CALL  RANDU(IX»IY»FY) 

IX=I  Y 

THETA-AL+ ( AU-AL ) *FY 

CALL  SSDOT  < 1 , 1 . ,  THETA , S » SDOT , NS » 0  > 

SUM1=0 « 

DO  60  1=1, NS 
60  SUM1=SUM1+S<I>**2 

SUM2-0. 

DO  70  1=1, NS 

70  SUM2=SUM2+S< I )#Y1 <  I ) 

AMP=SUM2/SUM1 

IF <AMP.LT . 1 .E-6)  AMP-1. E-6 
XJ=0. 

DO  80  1=1, NS 

80  XJ=XJ+<Y1<I)-S<I)*AMP>**2 

IF < X J . GE . XJO )  GO  TO  90 
IFLAG-1 
X ( NT ) =AMP 
X(NT*2)=THETA 
XJO-XJ 

90  CONTINUE 

IF < IFLAG.EO. 1 )  GO  TO  100 

X(NT)=l.E-6 

X<NT*2)=0. 

100  RETURN 
END 


C*  SUBROUTINE  SSDOT 

C 

C  COMPUTE  S ( T )  AND  ITS  DERIVATIVE  AT  NS  INSTANTS 

C 

C 

SUBROUTINE  SSDOT ( NT , AMP , THETA , S, SDOT , NS , IFLAG ) 

REAL*4  S  < 1 ) , SDOT < 1 > , AMP  < I > , THETA  <1 ) 

COMMON  /SAMPLE/BT , DT 

COMMON  /SPLINE/  TBEGIN, TDELTA,NPOINT,ARO< 1024 ), VAL1 < 1024 > , 
1  VAL2 < 1024  >  »SPLIN1 <1024,3), SPLIN2< 1024 , 3 ) 
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c 

C  IF  IFLAG-O*  SOOT  IS  NOT  COMPUTED 
C 

IF<NT.EQ.O>  00  TO  20 
N  -  NT42 
DO  10  1=1 »NS 
S(I)  -0. 

II  -  <1-1 >  4  N 
TI=BT+< 1-1 )4DT 
DO  10  K=1»NT 
T  =TI  -  THETAOO 

KT  =  IFIX<  <T-TBEQIN)/TDELTA)  +  1 
D  »  T-ARG(KT) 

51  »< <SPLIN1<KT.3>  *  D  +  SPLIN1<KT.2> > 

1  4D  +  SPLIN1  <KTt  1  > )  *  D  +  VALKKT) 

S<I>  =  S<I)  +  S14AMP<K> 

IF < IFLAG  .EQ.  0>  GO  TO  10 
K1  =  I1  +  K 
C 

C  SD0T<K1)  IS  THE  DERIVATIVE  OF  S  U.R.T.  B  INSTEAD  OF  A.  A-B442 
C 

SDOT <K1 )  *  SI  4  2.4S0RT<AMP(K>) 

K2  =  Kl+NT 

52  ■  < <SPLIN2(KTf3)4D  +  SPLIN2<KTr?> ) 

1  4D  +  SPLIN2<KT *  1 ) )  4D+  VAL2(KT> 

SDOT (K2)  — AMP(K)4S2 
10  CONTINUE 

RETURN 

20  NM=NT*24NS 

DO  30  1=1 'NS 
30  S(I>  =0. 

DO  40  1=1 »NM 
40  SDOT  < I ) =  0. 

RETURN 

END 


C4  SUBROUTINE  GMINV 

C 

C  MATRIX  INVERSION  ROUTINE 

C 

C 

SUBROUTINE  GMINV<NRrNC»A»U»MR»MT> 

C 

C  MATRIX  INVERSION  ROUTINE 

C  INPUT  ♦ 

c  NR.NC  =ROW  AND  COLUMN  DIMEN.  OF  A 

C  A  -MATRIX  TO  BE  INVERTED 

C  MT  PRINT  CONTROL  VARIABLE 

C  OUTPUT  : 

C  U  =  GENERALIZED  INVERSE  OF  A 

C  MR  =  RANK  OF  U 

C 

EXTERNAL  DOT 

DIMENSION  A( 1 ) »U< 1 ) .S(30> 

COMMON/MAI Nl/NDIM 
COMMON/ 1 NOUT /K I N  *  KOUT 
NDIM1  =  NDIM+1 
TOL  -  l.E-14 

ADV  -  l.E-24 

MR  »  NC 
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NRM1  -  NR-1 

TOLl  -  0. 

JJ  -  1 

DO  10  J-lrNC 

S<J>  -  DOT<NRiA< JJ)  »A< JJ) > 

IF  <S<J)  .OT.  TOLl)  T0L1-9<J> 

10  JJ  -JJ  +NDIM 

TOLl  -  ADV*  TOLl 
ADO  -  TOLl 
JJ-1 

DO  100  J-lfNC 
FAC  -  S(J) 

JH1  -J-l 
JRM  -  JJ+NRM1 
JCM  -  JJ+JM1 
DO  20  I-JJ» JCM 
20  U<I>  -  0. 

U<  JCM) -1 • 0 

IF (  J  .EQ.  D  GO  TO  34 
KK-1 

DO  30  K-lrJMl 

IF  <S( K )  .  EQ.  1.0)  GO  TO  30 
TEMP-  -DOT <  NR  *  A  <  J J ) » A  <  KK ) ) 

CALL  VADD<K»TEMP»U< JJ) f U<KK>  > 

30  KK-KK+NDIM 

DO  50  L-l »2 
KK-1 

DO  50  K-1.JM1 

IF  (S < K )  .EO.  0.)  GO  TO  50 
TEMP  — DOT<NR»A<JJ)fA(KK>> 

CALL  VADD <  NR . TEMP »A<JJ)»A(KK>) 

CALL  OADD<K » TEMP» U( JJ ) »U(KK) ) 

50  KK  =  KK+NDIM 

TOLl  -TOL  *FAC+ADV 

FAC  -  DOT<NR.A( JJ) »A< JJ) ) 

54  IF  < FAC  .OT.  TOLl)  GOT  0  70 

DO  55  I-JJ* JRM 

55  A(I)  -0. 

S<J>  =0. 

KK  =>  1 

DO  65  K-1.JM1 

IF  < S(K )  .EQ.  0.)  GO  TO  65 
TEMP  — DOT<K*U<KK) *U< JJ) > 

CALL  VADD  < NR* TEMP* A< JJ) *ACKK) ) 

65  KK  »  KK+NDIM 

FAC  -DOT< J*U< JJ) *U< JJ) ) 

MR  -MR-1 
GO  TO  75 
70  S<J>  -1.0 

KK-1 

DO  72  K-1.JM1 

IF  <8(K)  .EQ.  1. >  GO  TO  72 
TEMP-  -DOT(NR*A<JJ)*A<KK> > 

CALL  0ADD<K*TEMP*U<JJ) *U<KK> > 

72  KK-KK+NDIM 

73  FAC  -1 , /SORT < FAC) 

DO  80  I-JJ* JRM 

00  A( I )  -  A< I )  *  FAC 

DO  85  I-JJ* JCM 
85  U< I )-U< I )*FAC 

100  JJ-JJ+NDIM 

IF  <MR  .EQ.  NR  .OR.  MR.EC-NL)  GO  TO  120 
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c 

c 

c 

c 

c 


IF  (NT  .NE.  0)WRITE<K0UT f 1 10>NR*NC»MR 
FORMAT < 13 » 1HX » 12? 8H  Hi  RANKfI2> 

NEND  -NC*NDIH 
JJ-1 

DO  135  J-l.NC 
DO  125  I-1»NR 
II-I-J 
S(I)-0. 

DO  125  KK-JJ»NEND»NDIM 
Sd)-Sd)+A<  II+KKXU<KK> 

II-J 

DO  130  I-1>NR 
Udl>-Sd> 

II-II+NDIM 

JJ-JJ+NDIM1 

RETURN 

END 


SUBROUTINE  VADD<N.C1 » A.B) 
INPUT  i 

N  -  ARRAY  DIMENSION 
Cl  -  SCALAR 
A  -  NX1  CECTOR 
B  -  NX1  VECTOR 
OUTPUT  1 

A  =  NX1  VECTOR  SUM 
DIMENSION  A<1)'B<1) 

DO  1  I»1»N 


1  A(I)  =>A<  I  )  +C1*B<  I ) 

RETURN 
END 


FUNCTION  DOT <NRf ArB) 

C  INPUT  J 

C  NR  «  ARRAY  DIMENSION 

C  A  *  NRX1  VECTOR 

C  B  -  NR  X!  VECTOR 

DIMENSION  A( 1 ) >B(1 ) 

DOT  -0. 

DO  1  I-lrNR 

1  DOT-DOT  +  Ad  )*Bd ) 
RETURN 
END 
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